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ABSTRACT 

We present the first rates of flares from M dwarf stars in both red optical and near infrared (NIR) 
filters. We have studied -50,000 M dwarfs from the SDSS Stripe 82 area, and 1,321 M dwarfs from 
the 2MASS Calibration Scan Point Source Working Database that overlap SDSS imaging fields. We 
assign photometric spectral types from M0 to M6 using (r — i) and (i — z) colors for every star in 
our sample. Stripe 82 stars each have 50-100 epochs of data, while 2MASS Calibration stars have 
—1900 epochs. From these data we estimate the observed rates and theoretical detection thresholds 
for flares in eight photometric bands as a function of spectral type. Optical flare rates are found to 
be in agreement with previous studies, while the frequency per hour of NIR flare detections is found 
to be more than two orders of magnitude lower. An excess of small amplitude flux increases in all 
bands exhibits a power-law distribution, which wc interpret as the result of flares below our detection 
thresholds. In order to investigate the recovery efficiency for flares in each filter, we extend a two- 
component flare model into the NIR. Quiescent M0-M6 spectral templates were used with the model 
to predict the photometric response of flares from u to K s . We determine that red optical filters are 
sensitive to flares with w-band amplitudes >2 mag, and NIR filters to flares with Au >4.5 mag. Our 
model predicts that M0 stars have the best color-contrast for J-band detections, but M4-M6 stars 
should show the highest rate of NIR flares with amplitudes of AJ ~0.01 mag. Characterizing flare 
rates and photometric variations at longer wavelengths is important for predicting the signatures of 
M dwarf variability in next-generation surveys, and we discuss their impact on surveys such as LSST. 
Subject headings: stars: low-mass, stars: flare 



1. INTRODUCTION 

The study of flares on M dwarfs has a rich history, 
spanning more than seven d ecades (e.g. van Maanen 
1940[ | Joy fc Humason| [l~949[ ) . The similarity of these 
events to resolved flares on th e Sun was recognized many 
years ago (e.g . Lovell 1969). The seminal paper by 



are required to extend our understanding of stellar flare 
e mission physics. 



West et al. 



Kunkel ( 1970 1 laid the groundwork for our present un 



derstandmg of the broad multi-wavelength properties of 
M dwarf flares. The basic two-component flare spec- 
tral model (H reco mbination + w arm blackbody contin- 



2008 ) found that the fraction of M dwarfs 
that show evidence of magnetic activity as defined by Ha 
emission decreases with height above the Galactic plane. 
This is thought to be an age effect, with the older, less 
active stars lying further from the plane due to dynam- 
ical h eating and radial migr ation within the disk ( |West| 



et al. 2008 Loebman et al. 2011). Flares may play an 



uum) presented in |Kunkel (19701 is useful for estimating 



the total energy a nd stellar surface extent of flares in 



the oi 


atical 


et al. 


2010) 



flares, such as line emission from H, He I, He II, and Ca 
II H&K, nor the evolution of the spectral en ergy distri 
bution (SEP) as the flare cools and decays (Fuhrmeis 
ter et al.][2008 ) . A partnership of multi-wavelength pho 



important role in the main sequence stellar evolution for 
low-mass stars, and correlating M dwarf flare rates with 
mass and metallicity promises to yield an effective con- 
straint on the ages of stellar populations. 

Models and observations of the initial mass function 
indicate that M dwarfs are the primary stellar compo - 
nent in the Galaxy by number ( Bochanski et al.] |2010[). 
Recent l arge scale surveys, such as the Sloan Digital Sky 



tometric a nd spectroscopic time-domain ob serving cam 



paigns (e.g. Byrne 1989 Hawley et al 



radiative hydrodynamical modeling 



Survey ( |York et aJL.||2000[ hereafter SPSS) have pho 
tometrically identified millions of M dwarfs ( jBochanski] 
et al.|20 10). These stars make excellent tracers of nearby 
20031, and de tailed galactic structure due to their high number densities 
Alfred et al ||2006k ( jBochanski et al.|2007a[ ), and future survey mi ssions, like 



the L arge Synoptic Survey Telescope (LSST; Ivezic et al. 
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2008 1 will detect M dwarfs to nearly the full extent of the 
Galaxy, using longer wavelength bandpasses and deep 
repeat imaging. LSST will map the sky approximately 
1000 times over 10 years using ugrizy bands, and is ex- 
pected to detect flares both as variability in previously 
quiescent stars, and as transient emission from previously 
undetected sources. This foreground "fog" of character- 
istically blue transient events will prove a significant con- 
taminant to any statistical measurements of variability 
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fr om cosmological source s URau et al.||2008|. 



Stripe 

2007) to look 
in the 



Kow alski et al | ( |2009[ hereatter K09) us ed the SPSS 
82 time domain photometric data (Sesar et al. 
from 



for flares 
and g-bands. 



50,000 M dwarf light 
From these sparse light 



curves m tne u 
curves they were able to recover a flare rate consistent 
with those foun d from dedicated photometric monitoring 
campaigns (e.g. Lacy et al.|1976 ). The rate of flares was 
found to decrease with peak hare luminosity , and increase 
with st ellar spectral type (redder colors). |Hilton et al. 
(2010b) found serendipitous flares in the time-resolved 
SDKS spectroscopic data from ^38,000 M dw arfs across 
the entire SPSS Dat a Release 6 footprint (Adelman- 
McCarthy et al. 2008). They confirmed the increase of 



the hare duty cycle with later spectral type, and further 
investigated the Galactic height dependence of flaring ac- 
tivity, finding that the stars that flared were even more 
closely confined to the Galactic plane than the stars with 
Ha activity. 

The few studies of M dwarf flares in the near infrared 
(NIR) that have been previously attempted have not 
found a consistent relationsh ip to the optical emission. 
Rodono fc Cutispoto| ( |1988[ ), for example, measured a 
marginal anti-correlation between optical U-band and 
infrared K-band photom etry of two flares. Similarly, 



Panagi & Andrews 
during 26 M dwar: 



( 1995 ) found no K-band enhancement 
flares that were observed simultane- 
ously in the U-band. Observations of flares o n the S un 

|Xu et 



a more promising avenue. 



al. 



in the NIR offer 

(2006) studied the luminosity of solar flares at 1.6/xm, 
the H opacity minimum, and found a correlation with 
hard X-ray and white-light optical emission for two X- 
class flares. Their data indicated that the NIR band 
reached peak luminosity after the X-ray, and was re- 
sponding to secondary heating, possibly from chromo- 
spheric back- warming. 

Modern M dwarf flare studies are expanding both in 
wavelength and temporal coverage to gain insight into 
stellar flare physics, as well as the contaminating im- 
pact these events have on det ecting other sources of 
variability. Hawley et al. (2003[) found a correlation be- 



tween flare ene rgies in the optical, near-UV, and X-ray. 



Schmidt (2012) have observed flaring M dwarfs with NIR 



spectroscopy to characterize the correlation between NIR 
emission lines, such as the H Paschen series, and blue op- 
tical broadband flux during flares. Dedicated high time- 
cadence white-ligh t flare studies from the Kepler mission 



(Koch et al. 2010) will allow an order-of-magnitude im- 



provement m the stud y of flare frequency distribu tions 



( Walk owicz et al. 2011 ) and their possible evolution. Tof- 



flemire et al. (2011 (nave undertaken the most compre 



hensive simultaneous optical and NIR photometric ob- 
serving campaign for flares to date to determine the ef- 
fects of flares on exoplanet detection around M dwarfs. 
These high cadence and high precision data place up- 
per limits on the broadband NIR emission from several 
medium u-band flares from M3 and M4 stars with ^3.9 
milli-mag photometric precision. 

In this paper we seek to extend the analysis of K09 
using the time domain data from the 2MASS Calibration 
Point Source Working Database (hereafter Cal-PSWDB) 
as well as the SDSS Stripe 82 data. The combination 
of these data allows us to statistically characterize for 
the first time the rate of broadband M dwarf flares in 
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Fig. 1, — The distribution of the Cal-PSWDB tiles across the sky, 
with a few notable fields labeled. LMC fields are shown in green, 
SMC in red. Orange stars denote Cal-PSWDB tiles with matches 
to DR7. The SDSS footprint is shown for reference in grey, with 
the Stripe 82 data used shown in blue. 

red optical and NIR wavelengths. Our data are outlined 
in §2. Variability statistics are discussed in §3. Flare 
rates as a function of both spectral type and observed 
passband are presented in §4. A spectroscopic model of a 
M dwarf flare is given in §5, and the predicted sensitivity 
to finding flares in each band in §6. Finally a concluding 
discussion is given in §7. 

2. TIME DOMAIN DATABASES 

Our data are drawn from the SDSS and 2MASS pho- 
tometric time domain databases. Every M dwarf in 
our study has single-ep och SDSS Data Release 7 (DR7 
Abazajian et al. 2009) photometry, used to determine 
photometric spectral types and distances. Here we de- 
scribe the selection criteria and spatial matching of our 
time domain samples. 

2.1. 2MASS Cal-PSWDB 

Our near infrared data come from the 2MASS Cal- 
PSWDB. Full details are given in the online Explanatory 
Supplement]^] and Cutri et al. (2006). In summary, the 
repeat imaging was obtained to provide consistent pho- 
tometric calibration for the 2MASS survey over the four 
years of its operation. There were 35 regular calibra- 
tion tiles spread across the sky, with each tile spanning 
approximately 8'5 in RA, and 1° in Dec. One of these 
35 tiles was visited hourly during survey operations by 
the 2MASS telescopes, with simultaneous imaging in the 
J, H, K s bands. This sequence produced between 562 
and 3692 epochs for these calibration fields. Five ad- 
ditional tiles were imaged in the LMC and SMC with 
fewer epochs. A total of 191 million point source mea- 
surements are contained in the database. A pilot study 
with the Cal-PSWDB by |Becker et~aL] p008l) revealed a 
2.6 day low-mass eclipsing binary system~|Plavcha n et al.| 
( |2008b[) used the Cal-PSWDB data to search for periodic 
variability of low mass stars, fin ding several eclipsing bi 



nary candidates. Additionally, Plavchan et al 



(2008a) 



6 |http : //mm. ipac . caltechTed u/2mass/releases/allsky/doc/ 
seca4_l .html 
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found a periodic signature in Cal-PSWDB from a pos- 
sible triple young stellar object in the p Ophiuchus star 
forming region. 

We entered the Cal-PSWDB point source measure- 
ments into a MySQL database to allow rapid matching 
of these data to other datasets. The 191 million mea- 
surements were clustered into individual objects using 
the s patial clustering algorithm, OPTICS (Ankerst et al. 
1999J, creating light curves for 113,030 unique objects. 
Of the 35 primary Cal-PSWDB tiles 16 were coinci- 
dent with the SPSS Da ta Release 7 photometric catalog 
( Abazajian et al7] |2009[ , hereafter DR7), and we found 
1,445 matched point source objects using a 275 radius. 
Figure [I] shows the positions of the 40 Cal-PSWDB tiles 
as compared to the SDSS photometric footprint. An 
optical-NIR (g — K s ,i — J) color-color diagram for the 
matched point source objects is shown in Figure [2] with 
1,963 galaxies and 9,369 stars. The SDSS OBJ.TYPE 
flag was used to distinguish between stars and galaxies. 
We removed the 113 objects with unknown object type 
(OBJ_TYPE=0 ). The distribution of stars in our sample 
is well fit by the Covey et al. (2007) fiducial stellar locus. 

We imposed magnitude limits of r < 23 and z < 23, 
and standard SDSS photometric flag cuts to ensure good 
SDSS DR7 photometry. M dwarfs were selected and 
spectral types assigned using (r — i) and (i — z) colors 
and the covariance matrix technique as outlined in K09. 
We found 4,860 stars with photometric spectral types 
later than M0, with very few stars in the M7-M9 range 
due to their faint apparent magnitudes. The resulting 
4,731 M0-M6 stars had a total of 5.75xl0 6 point source 
measurements in Cal-PSWDB. Only 12 of ou r stars were 
in th e SDSS spectroscopic M dwarf catalog ( |West et al. 
2011 1 so we rely solely on the SDSS photometry for clas- 
sification. 

We restricted our analysis to epochs with the Cal- 
PSWDB flag PELQUAL = A, which required every 
JHK S measurement to have a signal-to-noise ratio of 
at least 10. To ensure that each DR7 source matched 
to only a single Cal-PSWDB object, we conservatively 
removed all objects with a spatial standard deviation in 
their Cal-PSWDB data of a a = ^ 



a 2 n ^ > 0'.'25. 



. T RA~ rw Dec 

This cut limited the effects of source confusion and de- 
blending, but also removed any high proper motion ob- 
jects from our sample. We additionally excluded epochs 
where photometric errors exceeded oj > 0.2, an > 0.25, 
or (Tk > 0.3. Our final Cal-PSWDB sample contained 
1,321 M dwarfs with ~2xl0 6 epochs, for an average of 
~1,900 good epochs per star. The selection of our Cal- 
PSWDB sample is detailed in Table [IJ 

2.2. SDSS Stripe 82 

In addition to the matched time domain JHK S data 
from Cal-PSWDB, we analyzed M dwarfs from the SDSS 
Stripe 82 repeat scan database (Sesar et al. 2007|) to 
obtain time domain data in the ugriz bands (Fukugita 
et al. |1996|). These data were previously used by K09 



to search for flares in the u- and <?-bands, and we adopt 
their sample as described in their Table 2. Briefly, these 
data contain 50,130 photometrically selected M dwarfs, 
with 1.9 million total photometric measurements. The 
stars were required to have good u-band magnitudes 
(it < 22), which limited the sample to relatively bright 
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Fig. 2. — A color-color diagram for the 10,000 point source ob- 
jects with matches in both DR7 and Cal-PSWDB. Stars, as de- 
termined by the SDSS OBJ.TYPE flag, are shown as red circles, 
galax i es as b lue. These two populations nicely separate. The |Covey| 
|et al.| l |2007[ l fiducial stellar locus is shown for reference, and follows 
closely tne stellar objects in our sample. 

M dwarfs. Only 429 of these objects have matches in the 
Cal-PSWDB database, due to the limited spatial over- 
lap. We measure the flare rate independently in each 
filter to build statistics using the entire population from 
each time- domain survey. 

The spatial location of the SDSS equatorial Stripe 82 
data is shown in Figure [T] in blue. Figure [3] gives an ex- 
ample of the eight-band light curves we have produced 
for a representative M0 star with both Stripe 82 and Cal- 
PSWDB photometry, illustrating the differences in time 
coverage and cadence between the two surveys. This star 
had 2,922 epochs of data in each of the JHK S bands, and 
70 epochs in each of ugriz. The time units for all of our 
SDSS and 2MASS data have been converted to barycen- 
tric dynamical time (BJDxdb)j & n d are presented as a 
modified barycentric Julian date (BJD-2400000.5). Fig- 
ure [4] shows the numbers of stars in each spectral type 
bin lor our Stripe 82 and Cal-PSWDB samples. 

3. THE * STATISTIC 

Discerning true variability from background noise re- 
quires sophisticated statistical techniques. Following 
K09, we utilize a modified Welch-Stetson variability in- 



son 


1993 


Stetson 


1996 



This statistic uses measure 



from the intrinsic scatter in a light curve. For the n- 
th measurement of a light curve in J- and iJ-bands, for 
example, $ is defined as: 



<f>jH(n) 



mj(n) - (mj) 



m H {n) 



{m H ) 



a H (n) 



(1) 



where m(n) is the apparent magnitude at epoch n, (m) 
the median magnitude over all epochs, and o~(n) the pho- 
tometric error at epoch n. When the magnitude increases 
or decreases for both bands simultaneously, $ yields a 
positive value. Conversely, when the magnitude from one 
band is higher than normal while the other is diminished, 
$ is negative. For purely random noise between the two 
filters, this index would yield a symmetrical distribution 
about $ = 0. Excesses of positive $ values represent cor- 
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TABLE 1 

Outline of sample selection cuts for the Cal-PSWDB data. 



Step 


Number 


Cut(s) 


Description 


1. 


11,445 (objects) 




Cal-PSWDB - DR7 Match 


2. 


5.75X10 6 (epochs) 4,731 (stars) 


<SpType < 6, OBJ.TYPE = 6 


M dwarfs only 


3. 


2.51x10 s (epochs) 1,321 (stars) 


PH.QUAL = A, 


High S/N Epochs 




1.14x10 s ( JH epochs) 


^A,Dec < 0"25 


Deblending 


4. 


J < 17, H < 16.5, a j < 0.2, a H < 0.25 


Good Photometry 


5. 


9.08x10 s (HK epochs) 


H < 16.5, K < 16, a H < 0.25, a K < 0.3 


Good Photometry 



2000 2001 2002 2003 2004 2005 



50000 



20.55 
18.07 
16.85 
16.32 
16.05 

1» 




51000 



J 
H 



52000 



14.88 
14.21 
14.08 



53000 



54000 



55000 



MJD T 



neously ( Skrutskie et al.||2006 ), allowing the photometry 
to be directly compared tor each 7.8 second exposure. 
The SDSS used a drift-scanning technique that produced 
separate 52.4 sec ond exposures in e ach filter taken in the 



order r,i,u,z,g (York et al. 20001. As a result the $ 



Fig. 3.— An eight-band light curve for SDSS J231820.76+010141.3, which had both Stripe 82 and Cal-PSWDB time-domain photometry. 
This object had a photometric spectral type of M0 assigned. Time has been converted for both surveys into BJD TDB. The median good 
magnitude for each filter is given beside each light curve. For reference the calendar year is given on the top-axis. 

related variability, from events such as flares or eclipses. 

We computed the $ statistic for all pairs of bandpasses 
"adjacent" in wavelength space (e.g. $„ fl , Qjh) for every 
epoch in our samples. Since the SDSS and 2MASS data 
were not obtained contemporaneously (see Figure [3]) we 
do not compute Q z j. Figure [5] shows the $ histograms 
for all objects in our Stripe 82 and Cal-PSWDB sam- 
ples. The black line indicates the $ > data, hereafter 
denoted $ + , while the blue line shows the |$ < 0| data, 
hereafter $ _ . 

As the sampling cadence of Stripe 82 and Cal-PSWDB 
surveys was much l onger than typ ical flare timescales 
of minutes to hours ( Moffett 1974 ) , our flare search re- 
lies on single-epoch outliers in the light curves. In or- 
der to act as an effective variability discriminant the 
two bandpasses making up the $ statistic must be im- 
aged nearly simultaneously to capture the same physical 
event. The 2MASS JHK S bands were imaged simulta- 



lO -1 



w 
w 

=fc 10 



1.0 r 



0.1 



Stripe 82 
Cal-PSWDB 



M0 Ml M2 M3 M4 M5 M6 M7 M8 M9 



values, for example, come from measurements which were 
separated by only ~1 minute, while the $> gr values are 
derived from data which were separated by almost 5 min- 
utes. We searched for flares in each of the six $ distribu- 
tions separately, but recognize the limitation that flares 
with timescales of only a few minutes may not have been 
imaged by all SDSS bands. As such, these rapid events 
would not yield large $ values, and may fall below our 
detection thresholds. 

To remove large $ values from the sample that resulted 
from a decrease in flux in both bands, for example due 
to photon loss from clouds or eclipses, we required all 
epochs in our investigation to have m < (m) for the 
bluer (shorter wavelength) filter in each of the $ + and 
$~ distributions. This removed approximately half the 
epochs in each distribution, as expected. As discussed in 
K09, the $~ distribution, also referred to as the "null" 
distribution, is representative of uncorrelated noise in the 
light curve. The 3>~ distribution was still well modeled 
by gaussian noise in all filters after the removal of half 
the epochs. Comparison of the <& + distribution with the 
$~ yields excesses in the $ + due to candidate flaring 
epochs. 

4. OPTICAL AND NIR, FLARE RATES 

In this section we explore the selection of flares from 
our light curves, and the luminosities of the flare candi- 
date epochs. 



Fig. 4. — The distribution of spectral types for our Stripe 82 and 
Cal-PSWDB samples. Due to the lack of stars in the M7-M9 bins 
we have only performed our analysis on M0 through M6 sub types. 



4.1. Observed Rates 

As discussed above, given our sparse sampling we as- 
sume that each flare is observed in only one epoch of our 
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Fig. 5. — Histograms of the two-band variability index, <J>, for each filter combination. The <J>+, where the flux in both filters is greater 
than their respective median values, are shown in black. The "I" - null distribution has been flipped about 3> = for comparison. We 
subtracted the null distribution from the positive distribution to measure the residual numbers of epochs with correlated increases in flux 
between bands. The <J> thresholds used in each distribution to define flare candidates are shown as vertical dashed lines. These are described 
in the text. 



data. A threshold value for $ + must be chosen from the 
null distributions to separate flare epochs from the qui- 
escent epochs with statistical fluctuations. If too high 
a threshold is chosen, we may lose valuable weak sig- 
nal in the redder bandpasses by disregarding real flare 
events. However, too low a threshold will add signifi- 
cant false positives to the resulti ng flare rates. Us ing 
the False Discovery Rate method ( Miller et al.||200lj ) on 



the $~ 9 distribution, K09 determined that a conserva- 
tive cutoff for flares was <E>+ ff > 100. This ensured less 
than 10% contamination from non-flaring epochs would 
be included in their final analysis. We performed the 
False Discovery Rate analysis to limit the contamina- 
tion to 10% in all six of our <I>~ distributions, and de- 
termined $ thresholds of (100, 100, 120, 500, 50, 50) 
for ($ ug , $ gr , $„, <$> lz , $ Jff , <5>hk s ) respectively. The $" 
distributions using 2MASS filters in Figure [5] decline 
rapidly at $ ~ 20, nearly an order of magnitude lower 
values than for the optical filters. This is indicative of 
the lower intrinsic stellar variability seen at these wave- 
lengths, and the larger photometric errors in the 2MASS 
survey. The $ thresholds are shown as dashed vertical 
lines in Figure [5] 

The number of flaring epochs for each filter combina- 
tion was very small compared to the total number of 
observations. Previous studies such as K09 have shown 
that flare rates for M dwarfs increase as a function of 
spectral type. To examine this trend we grouped our 
sample into three spectral type bins (M0-M1, M2-M3, 
M4-M6). The ratio of flare candidate epochs to the to- 
tal number of epochs for each bin is shown in Figure 

tAs expected, all the optical SDSS filter combinations 
ow an increasing numbers of flares with later spectral 
type (K09). There were considerably fewer flare epochs 
detected in the NIR passbands, although the and 
<&hk flare rates fall within the range of the SDSS data 
in Figure [6] This is likely due to the weak signatures of 
flares in these NIR bands, as discussed in §5. Table [2] 
presents the total numbers of flares in each band for the 



TABLE 2 

Summary of flare candidates in each spectral 
type bin per wavelength pair. 





Spectral 


Candidate 


Number 


Flare 






Type 


Flares 


Epochs 


Fraction 


ug 


MO-MI 


81 


1.5xl0 6 


5.4x10" 


5 


ug 


M2-M3 


116 


6.8xl0 5 


1.7x10" 


4 


ug 


M4-M6 


262 


l.lxlO 5 


2.3x10" 


3 


gr 


M0-M1 


722 


1.6xl0 6 


4.6x10" 


4 


gr 


M2-M3 


322 


7.2xl0 5 


4.5x10" 


4 


gr 


M4-M6 


251 


1.2xl0 5 


2.1x10" 


3 


ri 


MO-MI 


567 


1.6xl0 6 


3.6x10" 


4 


ri 


M2-M3 


283 


7.3 x10 s 


3.9x10" 


4 


ri 


M4-M6 


149 


1.2x10 s 


1.2x10" 


3 


iz 


MO-MI 


125 


1.6xl0 6 


8.0x10" 


5 


iz 


M2-M3 


100 


7.2 x10 s 


1.4x10" 


4 


iz 


M4-M6 


53 


1.2x10 s 


4.3x10" 


4 


JH 


MO-MI 


41 


6.3x10 s 


6.5x10" 


5 


JH 


M2-M3 


10 


3.0x10 s 


3.2x10" 


5 


JH 


M4-M6 


48 


2.0x10 s 


2.4x10" 


4 


HK 


M0-M1 


99 


5.1x10 s 


1.9x10" 


4 


HK 


M2-M3 


13 


2.4x10 s 


5.5x10" 


5 


HK 


M4-M6 


69 


1.7xl0 5 


4.1x10" 


4 



three spectral type bins. 

4.2. Flare Luminosities 

To compare our results in many photometric bands to 
those of previous studies, we must transform the flares 
from apparent magnitude enhancements to luminosities. 
The apparent magnitudes for every epoch and filter in 
our time-domain Stripe 82 and Cal-PSWDB samples 
were first converted to observed flux, in units of erg 



cm fi m , using the ca librations from Ivezic et al. 



and Cohen et al. 



(120031) respectively. To estimate 
the quiescent luminosity for each epoch, we obtained dis- 
tance estimates for our st ars from the {jj, r ,r — z) photo- 



metric parallax relation of Bochanski et al. (20101, using 



the (r — z) color from DR7. Errors in the computed lu- 
minosities are dominated by the distance uncertainties. 
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Fig. 6. — The fraction of flaring epochs, defined as the number 
of good epochs above the <E> threshold divided by the total number 
of good epochs, for each filter combination and spectral type bin. 



Many previous studies have shown that the frequency 
of fl ares decreases as a function of increasing flare energy 
Lacy et al.|1976[) or flare luminosity (K09) 
pass we selected the epochs identified as 



(e-g 



For each 
flare can- 



bandpass we selected the epochs 
didates in Table [l] and found a cumulative distribution 
in flare luminosity. The frequency per hour that these 
flares are observed was calculated by dividing the cumu- 
lative number distribution by the total observing time 
in hours for each filter, giving the flare frequency distri- 
butions plotted in Figure [7] These show for each filter, 
/, the number of flares pernour of observation, i/j, as a 
function of observed lumin osity, Lf . The YZ CMi flare 



frequency distribution from |Lacy et al. 



( 1976 ) is shown, 
and corresponds well with our optical distributions. 

The optical (ugriz) distributions are very similar, indi- 
cating that flares produce smooth continuous variations 
in these bandpasses. As with Figure |6l the JHK s -ba,nd 
flare frequency distributions in Figure]?] are highly under- 
sampled. NIR flares were recovered more than two orders 
of magnitude less frequently than for SDSS optical colors 
for the M4-M6 bin. 

The turnover at low luminosities seen in the SDSS fil- 
ters is due to incompleteness in the low-luminosity flare 
census. To explore whether the observed power-law dis- 
tributions (from log Lf— 29 to 31 erg s _1 ) extend to even 
lower luminosity flares than our $ thresholds recover, we 
moved the <I )+ detection thresholds to higher and lower 
values in the u-band. In Figure [lO] the power-law slope 
appears to extend to much lower luminosities as a result 
of the lower $ thresholds. We recovered many additional 
flare candidate epochs, but with a much higher rate of 
false positives. The high-luminosity end of Figure [10] 
(above log L u ^29.1 erg s _1 ) is insensitive to the choice 
of $ threshold. We present our flare frequency distribu- 
tions as observed, but note that it is possible to correct 
the frequencies for the efficiency in recovering flares as a 
function of luminosity (E. J. Hilton 2011 PhD Thesis). 
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Fig. 7. — The cumulative flare frequency distributions as a func- 
tion of observed luminosity in each passband for a sample spectral 
type bin. Errors in luminosity are dominated by photometric er- 
rors and assumptions in the distance. Incompleteness in recovering 
flares at the low-luminosity end results in the flattening seen clearly 
in each of the SDSS bands. 

4.3. Residual $ Distributions 

The $ distributions in Figure [5] indicate that there are 
many more epochs in the <& + distribution than in the $ _ , 
showing that the two bandpasses which make up each 
$ distribution preferentially increased in flux together. 
The excess between the <E> + and $~ is also evident at 
values of $ well below the flare candidate thresholds. By 
subtracting the <3> + from the 4>~ distributions for each 
of the six filter combinations in Figure [5] we are left with 
residual $ + distributions, as shown in Figure [8] for <& r j. 
These residual distributions follow a roughly power-law 
profile below the $ + thresholds, with a turn-over at low 
$ + amplitudes. We fit the residual distributions using 
the equation 



log 10 F = a + (p log 10 X) (10 



-tKx+s) 



(2) 



where Y is the $ + — $~ residual, a is the peak am- 
plitude for the residual, /3 a power-law slope at larger 
luminosities, 5 is the turnover $ between the power-law 
and the flattening at low $, 7 the exponential plateau 
slope at low and X is the $ + value. Figure [8] also 
illustrates where in the residual distribution each of the 
fit coefficients constrains. This functional form was cho- 
sen to simultaneously fit the residual power-law tail at 
large and the turnover at low $ values due to random 
noise dominating. The power-law coefficient f3 is plot- 
ted for each of the distributions in Fig 19} The /3 values 
ranged between -1.9 and -3.1. The higher slope in (5 in 
the $> gr fit may be due to the better photometric pre- 
cision in the g and r bands, compared to the u band, 
allowing those filters to be more sensitive to flares be- 
low the $ threshold. The Cal-PSWDB $hk s conversely 
shows the smallest observed residual, again suggesting 
that these longer wavelength filters are less affected by 
intrinsic stellar variability. However, all the fits are con- 
sistent with a value of (3 ~ —2 within the errors. 
The physical cause of this power-law slope, /3, at $ val- 
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Fig. 8. — The <I> r j_residual (black), and the corresponding fit 
(red) defined in Eqn[2] 

ues less than our flare cutoff threshold, may be the exten- 
sion of the increasing fl are frequency at lower luminosi- 
ties. Lacy et al. (1976) showed that the flare frequency 
distribution has a power- law shape, with a turnover at 
lower energies due to confusion with the quiescent vari- 
ability, and a fall-off at higher energies due to incom- 
pleteness. Their flare detection extended to lower energy 
flares than our analysis was sensitive to. This trend has 
bee n verified for M dwa rf flares more recently by K09 
and jHilton et al.| ( |2010b |. 

There is not a direct transformation to compare be- 
tween our observed $ residual distributions and the flare 



frequency distributions calculated by Lacy et al. ( 1976 1 
and others. Our residual $ distributions contain all the 
M sub-ty pes combined , for sta rs across a large range of 
distances. Hilton et al. (2010a) have shown that the flare 
frequency distribution slopes are similar, though not con- 
stant, for different spectral type bins, and between active 
and inactive stars 



Robinson et al. ( 1995 ) and Robinson et al. ( 1999 1 found 
that "microflares, very frequent low-energy rapid flares, 
were detectable using HST for two very active M dwarfs. 
By comparing the distribution of fluxes about the me- 
dian quiescent value in their light curves, the microflares, 
which are of too low amplitude to be detected individ- 
ually, could be detected statistically. They found that 
these microflare events followed a power-law distribution 
which was steeper than that of the larger resolved flares 
previously detected. Thus there are even greater num- 
bers of microflares than would be predicted by simply 
extending the power-law flare frequency to lower ener- 
gies. 

Our method of finding the residual $ distribution 
yields a similar result, by determining an asymmetry in 
the light curves via the <& statistic due to excess flux en- 
hancements. However we caution that we are not detect- 
ing the extremely low- amplitude microflares s een in the 
HST observations from |Robinson et al. (1995). Instead, 
the residual power in the <t> distributions are likely a re- 
sult of the underlying cacophony of low amplitude flares 
that are below our single-epoch detection limits in each 
bandpass, but which would be visible in high-precision 
continuous monitoring programs. It is especially notable 
that the significant asymmetric distributions were seen 
in all wavelengths studied (u-K s ), although we are sam- 
pling very different contrasts with the underlying quies- 
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Fig. 9. — The slope of the residual $ distributions, as defined in 
Eqnpl for each filter combination. The steeper slope at increasing 
wavelength is due to the lower intrinsic rate of variability within 
the photometric errors. 

cent stellar spectrum, which supports our interpretation 
of these events as flares. While the flare frequency dis- 
tributions in Figure [7] are under sampled for the NIR 
filters, the residual epochs with excess flux seen in the 
$ jh and &hks in Figure [5j and their power-law slopes 
in Figure [9j provide a statistically significant detection 
of small-amplitude flux enhancements. 

5. OPTICAL AND NIR FLARE MODEL 

We have employed a semi-anal ytical two-com ponent 
flare model, as first advocated by Kunkel (19701, to ex- 
plore the completeness of our search tor flares m every 
filter and to compare the o bserved flare luminosit ies be- 
tween filters. We extend the Kowalski et al. ( 2010 ) model 
which was used in blue optical wavelengths to red optical 
and NIR bandpasses. 

In order to simulate flux enhancements due to flares in 
each filter using the two-component model, we built qui- 
escent template M dwarf spectra in the optical and NIR. 
We utili zed flux calibrated spectra from the IRT F stellar 



catalog ( Cushing et aLp005[|Rayner et al.|2009[ ) for stars 
with spectral type M0 t hrough M6. The MO th rough 



M6 SDSS templates from Bochanski et al. (2007b) were 
used for the most of the optical regime. The 1K1F and 
Bochanski spectra were normalized and joined at 8500A. 
Since the Bochanski optical templates only go as short 
as 3825A, to fully cover the M dwarf spectr a in the it- 
band we used templates from Pickles (1998). As these 
are lower resolution than the Bochanski templates, we 
used them only from 2900A to 5000A. The Pickles and 
Bochanski templates are both published in normalized 
flux units, while the flux-calibrated IRTF data we used 
in the NIR are in observed flux. We found dist ances to 
each o f the IRTF M dwar f stand a rd stars from I Jenkinsl 



( |1952|), [van Altena et ah] ( |1995| , |van Leeuwenr f 2007 ) 



and" Jenkins et al.| ( |2009|) using the S lMBADl Iataba se, 
and standard radii from |N. Reid fc S. L. Hawley ( 2000 1 to 
convert the observed flux into surface luminosity for each 
star in the NIR. The normalized optical templates were 
then scaled to match at 8500A, and or resultant opti- 
cal/NIR spectral energy distributions for the M0 through 
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Fig. 10. — The cumulative flare frequency distribution, as in 
Figure [7] for the M-band using four different $ thresholds to select 
flare candidate epochs. Lower thresholds allow more contaminating 
epochs into our flare distribution, but the power-law slope seen at 
luminosities above ~ 10 29 erg s — 1 remains the same. 

M6 sub types are shown in Figure [Tl] 

The lower panel of Figure |T2] presents our M3 spectrum 
in arbitrary flux units i n black, with the ugr izJHK. 



transmission cur ves from |Fukugita et al7| ( 1996 ) and Co 
hen et al. (2003) shown for reference. The overall flare 



SFD shape was defined by a 10,000K blackbody con- 
tinuum. Previo us studies have found temperatures as 
low as ~8500K ffiajvley fc Fisher||1992| ) for flares, while 
Kowalski et al.| (]20 10D found evidence for te mperatures 



as high as 13,000Kr |Fuhrmeister et al.| ( |2008[) fit a black- 
body temperature of 20,000 K during the brief impul- 
sive phase of an enormous flare on CN Leo. These tem- 
peratures rapidly declined, however, with values around 
7,000-10,000 K for the remainder of the flare event. 

Shortward of 3646A, we added a second component of 
H Balme r continuum l i ne em ission, which was computed 
from the Allred et al. ( 2006 ) radiative hydrodynamic M 
dwarf atmosphere modeT The Balmer continuum was 
scaled to have 10 times larger surface area coverage than 
the blackbody component. The 10:1 scaling was approxi- 
mately an average of t hose used by previous studies, with 



Kowalski et al. (2010) using values from 3:1 to 16:1. Our 
idealized model does not include many detailed emission 
mechanisms that would alter these fluxes, such as line 
emission from other elements, or other levels of hydro- 
gen. Using measured s pectral line emis sion from M dwarf 
flares in the NIR from Schmidt (20121, we estimate that 
line emission would contribute an order of magnitude 
less than our predicted blackbody continuum emission 
in JHK S bands, and is therefore negligible for our anal- 
ysis. Continuum emission from higher levels of hydrogen 
in the flare atmosphere, and effec ts from chromo spheric 



([2006J, likely 



back- warming such as discussed in Xu et al. 
provide the greatest changes to the red optical and NIR 
flare SED. 

To investigate the effect of the flare model, we scale 
it by a surface coverage fraction, ranging from 5xl0 -7 
to 0.5, and add it to the M0-M6 template spectra. The 
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Fig. 11. — The combined 2900A to 2.5/jm spectral energy distri- 
butions for our M0-M6 templates. The ranges in wavelength over 
which each of the template sources contributes are shown at top. 

photometric response to flares in each filter was then cal- 
culated by conv olving the SPSS and 2MASS filter trans- 
mission curves ( |Fukugita et al.||1996| [Cohen et al.||2003| ) 
to the M0-M6 model spectra with and without the flare. 
In Figure 13 we show the transformations between the 
Aw-band and the A(grizJHK s ) band amplitudes for all 
the flare models we ran. 

The top panel of Figure [12] shows the changes in each 
photometric bandpass as computed by our model for 
a flare on an M3 star with 1% surface coverage for 
the Balmer continuum and 0.1% coverage for the black- 
body component. The resulting increases in flux for the 
ugrizJHK s bands were approximately (2.7, 0.62, 0.19, 
0.06, 0.03, 0.01, 0.01, 0.01) magnitudes respectively. 

The extremely blue response of the flare is due to a 
contrast effect between the intrinsically blue flare model 
and the underlying red stellar spectrum. Since the M3 
spectrum is much redder than the flare model, even large 
amplitude flares do not give a strong singal in the red- 
dest filters. However, K09 recovered a very bright flare 
in the Stripe 82 sample, visible in all five passband light 
curves, on an active M6 star. This flare had magnitude 
excursions of (5.50, 3.07, 2.34, 0.77, 0.29) for the ugriz 
bands respectively. Our M6 model, scaled to match the 
u-band excursion, predicted ugriz magnitude excursions 
of (5.5, 2.5, 1.2, 0.3, 0.1). In the 2MASS JHK S filters we 
predicted excursions of approximately (0.03, 0.01, 0.01). 
The coverage for this simulated flare was 0.06% of the 
stellar surface. As the two-component model does not 
include detailed effects, as discussed above, the model 
excursions agree with the observations to within a fac- 
tor of two in flux. The K09 flare may also have evolved 
in temperature over the co urse of the SDSS observations 
(Fuhrmeister et al. 2008). We allowed the blackbody 
temperature within our model to vary from 8,000K to 
13,000K. Lower temperatures yielded only a marginally 
improved fit for the r-band, and required a larger effec- 
tive surface area for the flare. 

6. MODEL FLARE RECOVERY 

Using our flare model with the M3 spectral template, 
we investigated the efficiency of recovering flares from 
our sample with the $ statistic. We first computed the 
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Fig. 12. — Bottom: Our two-component flare model (dashed line) shown against the quiescent M3 optical-NIR stellar SED (black line). 
The ugrizJHK s filter passbands are shown for reference. Top: The resulting change in magnitude for each filter from a sample flare with 
1% surface coverage on an M3. 



sigma-clipped mean photometric error profiles as a func- 
tion of observed magnitude for each of the ugrizJHK s 
bands in increments of 0.1 mag, shown in Figure [14] For 
each of the M3 stars in our NIR time domain sample, 
we simulated 200 artificial flares using our model from 
§5, with a logarithmically increasing surface area cover- 
age fraction from 10 -7 to 10 _1 . The mean error profiles 
were then used to select an appropriate photometric error 
for the resultant apparent magnitudes of each simulated 
flare in all eight passbands. Finally, we calculated the $ 
values for each simulated flare on every star. 

The result was a "$ profile" , which was a function of 
both the quiescent apparent magnitude of the star, and 
the flare amplitude. These represent the two competing 
physical quantities, the distance to the star and the size 
of the flare, which conspire to limit our detection of small 
flares or flares on distant stars. This envelope is shown in 



Figure 15 as a function of quiescent magnitude on the x- 
axis, with $ contour values increasing from black to red. 



The $ profile was smoothed with a kernel of 0.1 magni- 
tudes. Because flares are typically characterized by their 
enhancement of blue flux, the corresponding predicted 
u-band excursion from our model is shown on the right 
axis in blue for the grizJHK s bands. 

Using our described <& thresholds from §4, we found 
that small flares in the ti-band (Am ~ 0.5mag) would be 
recovered for every star in the sample with a quiescent 
magnitude of u < 20.5. The completeness is a direct 
function of the photometric error profil e, w hich parallels 
the $ threshold curves shown in Figure 15 Interestingly, 
note that in the griz bands in Figure 15 the $ profiles 
"turn up" for bright stars, since the error profiles be- 
gin to increase as stars approach the saturation limit for 
the SDSS. Consider a star with a moderate quiescent 
luminosity in the u-band (e .g. u ~ 17). A very large 



flare such those observed i n Hawley & Pettersen ( 1991 1 
and Kowalski et al. ( 2010[ ) would result in an observed 
u-band magnitude too bright for the SDSS photometry, 
and would therefore be excluded from our analysis. This 
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Fig. 13. — The transformation between the predicted u-band 
response and the grizJHK s -band response for each spectral type 
template. 
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Fig. 14. — Mean photometric error as a function of apparent 
magnitude for all the Stripe 82 and Cal-PSWDB measurements in 
our database. 
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Fig. 15. — Lines of constant $ as a function of both mean apparent magnitude (distance) and change in ma gnitude (flare amplitude) 
in each filter, using the M3 spectral flare model above, and the mean photometric error profiles from Figure |T4| The predicted change in 
■tt-band which corresponds to the change in each other passband is shown on the right axis in blue. Contours from black to red are <E>=(50, 
100, 300, 700, 1200, 2000, 5000, 8000). These represent our detection thresholds for a given flare amplitude as a function of apparent 
magnitude in these surveys. The dotted-line threshold contours for 2MASS filters are <J>=5 and 10. 



would be possible in the g, and even r-band, for bright 
stars with very large flares. Incompleteness is therefore 
also expected at the brightest flare luminosities, although 
these events are rare with duty c ycles estimated to be ~1 
per month on most active stars ( Hilton et al.|2010b ). 

The erro r profiles in the <?, r, and i Alters plotted in 
Figure 14 are best suited to monitor small amplitude 
var iati ons such as eclipses and spot modulation. Figures 
[15] and [13] clearly show that we are only completely sam- 
pling flares larger than Au ~ 2 in these bands for most 
of our sample, and thus the u-band is the most suitable 
to search for the frequent low-amplitude flares in these 
stars. 

As expected from the model discussion above, only 
very large flares would be detectable in the NIR pass- 
bands for the brightest stars in our sample using the 
2MASS photometric error profile. Our model predicts 
that a flare with a nearly ~ 5 magnitude increase in the 
u-band would be required to detect any significant in- 
crease in flux in the Cal-PSWDB for a nearby M3 type 
star. No simulated flares on M3 stars showed enough 
increase in flux to register as a $ = 50 response in the 
JHK S bands, however. The flare emission in the NIR 
must therefore be greatly influenced by physics not in- 
cluded in our two-component model, and definitive si- 
multaneous detection of flares in the NIR and optical is 
needed to calibrate this model. 

6.1. Optimal NIR Flare Detection 

Our analysis has extrapolated a well-known two- 
component flare model into a new wavelength regime, 
and we therefore provide a prediction for the "best-case" 
to detect M dwarf flares in the NIR. An enhancement 
of AJ = 0.01 magnitudes is easily observable for many 
nearby M dwarfs with present day ground based differ- 
ential photometry. Using the Ait to A J flare amplitude 
transformations shown in Figure [13} we estimated the Au 
enhancement for each spectral type required to achieve 



AJ = 0.01. For spectral types M-M6 this yielded a 
Au =(1.4, 1.8, 2.0, 2.5, 3.3, 4.0, 4.5) mag. The changing 
shape in the underlying stellar SED in the NIR indicates 
that J-band flares are detectable from smaller Am flares 
at earlier spectral types. 

The frequency of flaring has been previously found to 
increase with spectral type, such as in our Figure [6] The 
characteristic lifetime of mag netic activity is a lso nighly 
correlated with spectral type ( |West et al. 2008 ) , and later 
spectral types M dwarfs are expected to show not only 
more, but larger amplitude flares (Hilton 2011 PhD The- 
sis). This highlights the contrasting effects at work in 
detecting broadband NIR flares: Persistently active M0 
stars that show large amplitude flares are exceedingly 
rare, as compared to those at later spectral types. The 
higher bolometric luminosity of earlier spectral type stars 
additionally means much larger energy flares are needed 
to achieve these large amplitudes. 

Figure 16 shows the observed M-band flare frequency 
distribution for each spectral type bin as a function of 
Au magnitude. The total rate of flaring increased with 
spectral type, and the slopes were very similar for all 
the distributions, in agreement with results from E. J. 
Hilton (2011 PhD Thesis). A linear fit was found for each 
logio^(Au) distribution between 0.4 < Au < 3.0 to avoid 
false-positives at the low-amplitude end, and incomplete- 
ness at high Au amplitudes. The predicted frequency 
at which we expect to observe a AJ = 0.01 mag flare, 
from the corresponding Au amplitudes, at each spectral 
type was then calculated from these fits, and is shown 
in Figure 17 The optimal spectral type to search for 



these NIR flare detections is in our latest type stars, M6. 
The slight up-turn at M0 is due to the more favorable 
color-contrast between the flare and the star. The over- 
all trend from M1-M6, however, shows that the strongest 
effect for all spectral types is the total flaring duty cy- 
cle. This provides one possible reason that the recovered 
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Fig. 16. — The observed flare frequency distribution for u-band 
flares at each spectral type. Each distribution was fit with a linear 
polynomial in the range 0.4 < An < 3.0. 

flare rates are higher than expected for MO-MI stars in 
Figure [6] A new method for calculating the duty cycle 
of a given amplitude excursion, known as the Magnitude 
Likelihood Distribution, has been recently put forward 
by E. J. Hilton (2011 PhD Thesis). Using these Mag- 
nitude Likelihood Distributions we predict a peak duty 
cycle of ~ 3 x 10~ 4 for the magnetically active M4 and 
M5 stars. For the 54 sec exposures of SDSS, this is equal 
to a frequency of v ps 3 x 10~ 2 hr _1 , which is consis- 
tent with our estimation. We would encourage future 
NIR flare observers to monitor bright /nearby targets in 
the latest spectral types possible, at least M4 or later. 
Promising stars might include GJ 1061, CN Leo (Wolf 
359), and especially the highly active YZ CMi. 

7. SUMMARY AND CONCLUSIONS 

We have shown that the signatures of flares from M 
dwarfs are not confined to the bluest bandpasses in large 
scale photometric surveys. There is significant response 
from medium to large flares in gri bands, and the largest 
flares are observable in z and in some cases into the NIR 
passbands. This is the first statistical characterization of 
M dwarf flares in the NIR regime. 

A self-consistent result between all eight filters shows 
that while the frequency of flares decreases with increas- 
ing luminosity in all filters, the total number of flares 
recovered decreases with wavelength. This suggests that 
like our simplified flare model the underlying SED for 
flares of all sizes is intrinsically blue. Furthermore, the 
total rate of flares in all optical bands appears to increase 
with increasing spectral type, but low-number statistics 
prevent us from robustly exploring this distribution at 
NIR wavelengths. The slope of the flare frequency dis- 
tribution is found to be insensitive to the flare detection 
threshold below ~ 10 29 erg s _1 . 

We have also presented an implementation of the two- 
component spectral flare model from |Kowalski et al.| 
(2010) to make predictions of the red-optical and NIR 
response of flares in M dwarfs. These conversions will 
be useful in future broadband photometric studies of M 
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Fig. 17. — The rate per hour at which a AJ = 0.01 mag am- 
plitude flare is expected to be seen on each spectral type, based 
on our model and measured flare rates. The numbers below each 
point are the Au amplitudes that our flare model predicts from a 
AJ = 0.01 mag flare at each spectral type. 

dwarfs to predict the contamination from flares across a 
wider wavelength space. Flares in the J-band are pre- 
dicted to be detected most frequently in spectral types 
from M4 to M6, with rates of - 2 x 10~ 2 hr" 1 . 

Medium sized flares larger than 0.5 magnitudes in u 
were found to occur at a rate of one per few thousand 
SDSS exposures for M dwarfs. Low-amplitude flare rates 
are limited by the photometric precision of the survey in 
question, and LSST will detect even smaller amplitude 
flares from these stars, and will increase the rate at which 
variability from these sources is detected. LSST will 
also employ a y-filter, with a central wavelength between 
that of z and J. Using the "y4" filter curve for LSST, 
our two-component flare model from Figure [12] predicts 
a flare amplitude of At/4 ~ 0.02 mag. The increased 
photometric sensitivity will therefore allow LSST to rou- 
tinely detect large flares in this NIR bandpass for bright 
stars, providing useful constraints on the SED structure 
of flares at long wavelengths 

Finally, we have shown that the intrinsic SED of a 
flare, convolved with the photometric characteristics of 
a given survey, can greatly affect the sample of flares 
recovered. The cadence and integration times will also 
affect the recovery of flares. Short exposures such as in 
2MASS capture the near instantaneous luminosity of an 
event, while longer exposure times may encapsulate en- 
tire flare events. We anticipate a day when this problem 
can be inverted, using orders of magnitude larger samples 
of flare candidate epochs in many bands, as well as spec- 
tra, from next-generation surveys to probe the underly- 
ing flare SED (or its deviations) across the population of 
M dwarfs in the Galaxy. 
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